// panda.cpp : Defines the entry point for the console application.
//

#if _MSC_VER > 1000
#pragma warning( disable :  4511 4512 4786)
// nasty hack for VC++ non-compliance with standard C++ for loop declaration scope
#define for if (0); else for
#endif

#include <stdio.h>
#include <math.h>
#include <string>
#include <set>
#include <iostream>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <iterator>
#include <sstream>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <boost/tuple/tuple.hpp>
#define _MY_STL_MACROS
#include "components/stlhelper.h"
#include "components/perlutil.h"
#include "components/popen/popen.h"
#include "components/progress_indicator.h"
#include "components/ultoa.h"
#include <assert.h>
#include "blosum_data.h"
#include "components/portable_timer.h"
#include "get_arg.h"

using std::string;
using std::stringstream;
using std::cerr;
using std::cin;
using std::vector;
using std::fill;
using std::cout;
using std::ifstream;
using std::ofstream;
using std::ostream;

#define BLOSUM80
#ifdef BLOSUM60
#define GAP_EXIST_PENALTY -11
#define BLOSUM_MATRIX blosum62
#else
#ifdef BLOSUM80
#define GAP_EXIST_PENALTY -10
#define BLOSUM_MATRIX blosum80
#endif
#endif
// global variables!!
// This must be a perl programme in disguise!!
namespace
{
	string my_path;
	string tmp_file_path;
	string yn00_exe_path;
	string paml_path;
	char lookup_table[256][256];
	char gap_penalty[256][256];
}

string GetPathPart(const string& path)
{
	string::size_type p1, p2, p;
	p1 = path.rfind('/');
	p2 = path.rfind('\\');

	// neither found : return this directory
	if (p2 == p1)
		return ".";

	if (p1 != string::npos)
	{
		if (p2 != string::npos)
			p = std::max(p1, p2);
		else
			p = p1;
	}
	else
		p = p2;
	return path.substr(0, p);
}




//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

//		t_pg_copy_guard

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

class t_pg_copy_guard
{
	t_pg_cursor&	pg_db;
	std::string		table_name;
	std::ostream&	CERR;
	bool			proceeding;
	bool			success;
public:
	std::string		buf;

	t_pg_copy_guard(t_pg_cursor& pg_db_, const std::string& table_name_, std::ostream& CERR_);
	~t_pg_copy_guard(){end_copy();}
	void cached_write();
	void write();
	operator bool() const { return success; }
	bool end_copy();

};


t_pg_copy_guard::t_pg_copy_guard(t_pg_cursor& pg_db_, const string& table_name_, ostream& CERR_)
	:	pg_db(pg_db_), table_name(table_name_), CERR(CERR_),
		proceeding(false), success(false)
{
	string copy_cmd = "COPY " + table_name + " FROM STDIN";
	if (pg_db.Exec(copy_cmd.c_str()) != PGRES_COPY_IN)
	{
		CERR	<< "\nERROR!! There was an error COPYING hash data to the postgreSQL "
				<< "table " << table_name << "\n"
				<< "[" << copy_cmd << "]\n"
				<<	pg_db.ErrorMessage()
				<<	std::endl;
		pg_db.PutLine ("\\.\n");
		pg_db.EndCopy();
	}
	else
		proceeding = success = true;
}
void t_pg_copy_guard::cached_write()
{
#ifdef NDEBUG
	if (buf.length() > 32000)
#else
/*	DEBUG
	std::cerr << "\n" << buf;
*/
//	if (buf.length() > 1)
#endif
	{
		write();
	}
}

void t_pg_copy_guard::write()
{
#ifndef NDEBUG
/*	DEBUG
	std::cerr << "\n" << buf;
*/
#endif

	if (buf.length())
	{
		pg_db.PutLine (buf.c_str());
		buf.clear();
	}
}
bool t_pg_copy_guard::end_copy()
{
	write();
	if (!proceeding)
		return false;
	proceeding = false;
	pg_db.PutLine ("\\.\n");
	if (pg_db.EndCopy())
	{
		CERR	<< "\nERROR!! There was an error COPYING hash data to the postgreSQL "
				<< "table " << table_name << "\n"
				<<	pg_db.ErrorMessage()
				<<	std::endl;
		return success = false;
	}
	return success = true;
}








//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
//	loci code
//
struct t_segment
{
	unsigned beg;
	unsigned end;
	t_segment(unsigned beg_, unsigned end_)
		: beg(beg_), end(end_){}
	t_segment(const string& data)
	{
		beg = atoi(data.c_str()) - 1;
		end = atoi(data.c_str() + data.find(':') + 1) - 1;
		if (beg  > end)
			beg = end = 0;
	}
};
typedef vector<t_segment> t_loci;



//________________________________________________________________________________________

//	merge_overlapping_loci

//________________________________________________________________________________________
void merge_overlapping_loci(t_loci& loci,
							t_loci& complement_loci,
							unsigned aln_beg, unsigned aln_end)
{
	vector<char> bits(aln_end, '\x0');
	typedef vector<char>::iterator iter_pos;

	// mark bits of all masked regions
	for (unsigned i = 0; i < loci.size(); ++i)
	{
		fill(	bits.begin() + std::min(loci[i].beg, aln_end),
				bits.begin() + std::min(loci[i].end + 1, aln_end),
				0x1);
	}

	// non-overlapping masked loci
	t_loci loci2;

	// get masked regions
	iter_pos first = find(bits.begin() + aln_beg, bits.end(), 0x1);
	while (first != bits.end())
	{
		iter_pos last = find(first+1, bits.end(), 0x0);
		loci2.push_back(t_segment(first - bits.begin(),
								   	   last  - bits.begin() - 1));
		first = find(last, bits.end(), 0x1);
	}
	loci2.swap(loci);

	// get non-masked regions
	first = find(bits.begin() + aln_beg, bits.end(), 0x0);
	while (first != bits.end())
	{
		iter_pos last = find(first+1, bits.end(), 0x1);
		complement_loci.push_back(t_segment(first - bits.begin(),
												 last  - bits.begin() - 1));
		first = find(last, bits.end(), 0x0);
	}
}

//________________________________________________________________________________________

//	subtract_overlapping_loci

//________________________________________________________________________________________
void subtract_loci(	const t_loci& original_loci, const t_loci& remove_loci,
					t_loci& resulting_loci,
			   		unsigned aln_beg, unsigned aln_end)
{
	resulting_loci.clear();
	vector<char> bits(aln_end + 1, '\x0');
	typedef vector<char>::iterator iter_pos;

	// mark bits of all loci
	for (unsigned i = 0; i < original_loci.size(); ++i)
	{
		fill(	bits.begin() + std::min(original_loci[i].beg, aln_end),
				bits.begin() + std::min(original_loci[i].end + 1, aln_end),
				0x1);
	}

	// remove bits of all loci
	for (unsigned i = 0; i < remove_loci.size(); ++i)
	{
		fill(	bits.begin() + std::min(remove_loci[i].beg, aln_end),
				bits.begin() + std::min(remove_loci[i].end + 1, aln_end),
				0x0);
	}

	// get resulting loci
	iter_pos first = find(bits.begin() + aln_beg, bits.end(), 0x1);
	while (first != bits.end())
	{
		iter_pos last = find(first+1, bits.end(), 0x0);
		resulting_loci.push_back(t_segment(first - bits.begin(),
			   						   last  - bits.begin() - 1));
		first = find(last, bits.end(), 0x1);
	}
}


//________________________________________________________________________________________

//	 extract_seq_pair_with_loci

//________________________________________________________________________________________
inline string extract_seq_with_loci(const std::string& seq, const t_loci& loci)
{
	string nrvo;
	for (unsigned i=0; i < loci.size(); ++i)
		nrvo += seq.substr(loci[i].beg,  loci[i].end + 1 - loci[i].beg);
	return nrvo;
}
inline void extract_seq_pair_with_loci(vector<string>& sequences,
										const std::string& seq1,
										const std::string& seq2,
										const t_loci& loci)
{
	// add sequence for mouse
	string seq_region;
	for (unsigned i=0; i < loci.size(); ++i)
		seq_region += seq1.substr(loci[i].beg,  loci[i].end + 1 - loci[i].beg);
	sequences.push_back(seq_region);

	// add sequence for human
	seq_region.clear();
	for (unsigned i=0; i < loci.size(); ++i)
		seq_region += seq2.substr(loci[i].beg,  loci[i].end + 1 - loci[i].beg);
	sequences.push_back(seq_region);
}

//________________________________________________________________________________________

//	 extract_cdna_pair_with_loci

//________________________________________________________________________________________
struct op_prot_loci_to_dna
{
	void operator() (t_segment& a)
	{
		a.beg *= 3;
		a.end = (a.end + 1) * 3 - 1;
	}
};

inline void extract_cdna_pair_with_loci(vector<string>& sequences,
										const std::string& seq1,
										const std::string& seq2,
										const t_loci& prot_loci)
{
	// prot loci to cdna loci
	t_loci loci = prot_loci;
	for_each (beg_to_end(loci), op_prot_loci_to_dna());

	// share processing with protein function
	extract_seq_pair_with_loci(sequences, seq1, seq2, loci);
}

//
//	loci code
//
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888






//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	bitscore

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//________________________________________________________________________________________

//	get_seq_bitscore

//		calculate blosum raw bitscores for one stretch of a sequence pair

//________________________________________________________________________________________
// helper struct
struct op_calc_bitscore
{
	char	(*gap_penalty)[256];
	char	(*lookup_table)[256];
	int	bitscore;
	char prev_a;
	char prev_b;
	op_calc_bitscore(char gap_penalty_[256][256], char lookup_table_[256][256]):
			gap_penalty(gap_penalty_), lookup_table(lookup_table_),
			bitscore(0), prev_a('-'), prev_b('-'){}
	void operator()(char a, char b)
	{
		// lookup bitscore
		bitscore += lookup_table[a][b];

		// gap opening penalty
		bitscore += gap_penalty[prev_a][a];
		bitscore += gap_penalty[prev_b][b];
		// save current char to calculate gap opening penalty
		prev_a = a;
		prev_b = b;

	}
};

int get_seq_bitscore(const std::string& seq1, const std::string& seq2,
							unsigned beg, unsigned end, unsigned& len)
{
	if (!seq1.length())
		return 0;
	// skip past initial '-'
	beg = seq1.find_first_not_of('-', beg);
	beg = seq2.find_first_not_of('-', beg);

	// skip past final '-'
	end = seq1.find_last_not_of('-',  end);
	end = seq2.find_last_not_of('-',  end);
	if (beg == string::npos ||
		end == string::npos || end < beg)
		return 0;
	assert(beg >= 0);
	assert(end < seq1.length());
	assert(seq2.length() == seq1.length());
	len += end + 1 - beg;
	return (for_each_pair(	seq1.begin() + beg, seq1.begin() + end + 1,
							seq2.begin() + beg,
							op_calc_bitscore(gap_penalty, lookup_table))
			).bitscore;
}

//________________________________________________________________________________________

//	get_all_regions_bitscores

//________________________________________________________________________________________
// helper function
void print_bitscore_over_loci(ostream& os, const std::string& seq1, const std::string& seq2,
								const t_loci& loci)
{
	double bitscore = 0.0;
	unsigned len = 0;
	for (unsigned i=0; i < loci.size(); ++i)
		bitscore += get_seq_bitscore(seq1, seq2, loci[i].beg, loci[i].end, len);
	if (len)
		os << "\t" << (bitscore / len);
	else
		os << "\t\\N";
}

void get_all_regions_bitscores(	ostream& os, const string& prot_m, const string& prot_h,
								const vector<t_loci>& vec_loci)

{
	// bitscore for each region
	for (unsigned i = 0; i < vec_loci.size(); ++i)
		print_bitscore_over_loci(os, prot_m, prot_h, vec_loci[i]);
}





//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	percentage_id

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


//________________________________________________________________________________________

//	get_all_regions_percentage_id

//________________________________________________________________________________________
//	helper functions
struct t_count_identical
{
	explicit t_count_identical(unsigned i = 0): count(i){}
	unsigned count;
	void operator() (char a, char b)
	{
		if (a == b)
			++count;
	}
};

void get_all_regions_percentage_id(ostream& os, vector<string>& proteins)
{
	// percentage id = -1.0 for empty strings
	for (unsigned i = 0; i < proteins.size(); i += 2)
	{
		if (proteins[i].length())
		{
			os << "\t" << (
							 for_each_pair(	beg_to_end(proteins[i]),
									proteins[i+1].begin(),
									t_count_identical()).count
									*100.0 / proteins[i].length());
		}
		else
			os << "\t\\N";
	}
}


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	KaKs

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//	helper function for below
// returns whether isnan / isinf /is 0.0
inline bool print_valid_double_pg_field(double num, ostream& os)
{
	if (isnan(num) || isinf(num) || num > 50)
	{
		os << "\t\\N";
		return false;
	}
	else
	{
		os << "\t" << num;
		return true;
	}
}

//________________________________________________________________________________________

//	print_sequence_ka_ks

//________________________________________________________________________________________
bool print_sequence_ka_ks(ostream& os, const string& cdna_m, const string& cdna_h,
							ofstream& err_strm)
{
	if (!cdna_m.length())
	{
		os << "\t\\N\t\\N\t\\N";
		return true;
	}
	// remove any existing results file before running analysis
	remove( "2YN.dN");
	remove( "2YN.dS");

	ofstream seq_strm ("seq.tmp");
	if (!seq_strm)
	{
		cerr << "\n\nError!!!\n\tCould not open /seq.tmp\n";
		return false;
	}

	seq_strm << "2 " << cdna_m.length() << "\n";
	seq_strm << "a  \n" << cdna_m << "\n";
	seq_strm << "b  \n" << cdna_h << "\n";
	seq_strm.close();

	{
		Popen Run_yn00;
		FILE *FH_YN00_ERR;
		FILE* FH_YN00 = Run_yn00.open_stderr(yn00_exe_path.c_str(), &FH_YN00_ERR, "r");

		if (!FH_YN00)
		{
			cerr << "\n\nError!!!\n\tCould not run " << yn00_exe_path <<"\n";
			return false;
		}

		// log errors
		char err_buf[1024];
		while (size_t bytes = fread(err_buf, 1, sizeof (err_buf) , FH_YN00_ERR))
			err_strm.write(err_buf,  bytes);
		Run_yn00.close_stderr(FH_YN00, FH_YN00_ERR);
	}

	// get results
	ifstream f_Ka ("2YN.dN");
	ifstream f_Ks ("2YN.dS");
	f_Ka.ignore(1000, 'b');
	double Ka;
	f_Ka >> Ka;
	f_Ks.ignore(1000, 'b');
	double Ks;
	f_Ks >> Ks;
	if (! f_Ka || !f_Ks)
	{
		os	<< "\t\\N\t\\N\t\\N";
		return true;
	}

	bool valid1 = print_valid_double_pg_field(Ka, os);
	bool valid2 = print_valid_double_pg_field(Ks, os);
	if (valid1 && valid2 && Ks != 0.0 && Ks != -0.0)
		print_valid_double_pg_field(Ka/Ks, os);
	else
		os << "\t\\N";
	return true;
}

//________________________________________________________________________________________

//	get_all_regions_ka_ks

//________________________________________________________________________________________
bool get_all_regions_ka_ks(ostream& os, ofstream& err_strm, const vector<string>& cDNAs)
{

	for (unsigned i = 0; i < cDNAs.size(); i += 2)
		if (!print_sequence_ka_ks(	os,
									cDNAs[i],
									cDNAs[i+1],
									err_strm))
			return false;

	return true;
}



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	matrix

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//________________________________________________________________________________________

//	get_all_regions_matrix

//________________________________________________________________________________________
struct op_calculate_matrix
{
	op_calculate_matrix(short freq_lookup_table_[256][256],	short* freq_matrix_):
			freq_lookup_table(freq_lookup_table_), freq_matrix(freq_matrix_){}
	short	(*freq_lookup_table)[256];
	short* 	freq_matrix;
	void operator() (char a, char b)
	{
		freq_matrix[freq_lookup_table[a][b]]++;
	}
};

void get_all_regions_matrix(	ostream& os,
								vector<string>& sequences,
							  	short freq_lookup_table[256][256])
{
	short freq_matrix[401];

	for (unsigned i = 0; i < sequences.size(); i += 2)
	{
		memset(freq_matrix, 0, sizeof (freq_matrix));
		for_each_pair(	beg_to_end(sequences[i]), sequences[i+1].begin(),
						op_calculate_matrix(freq_lookup_table, freq_matrix));

		//	only print out differences in the two sequences
		os << "\t";
		for (unsigned j = 0; j < 400; ++j)
			if (freq_matrix[j] != 0 && (j / 20 != j % 20))
				os << valid_aa_chars[j / 20] << valid_aa_chars[j % 20]
					 << freq_matrix[j] << ",";

		//	only print out identical
		os << "\t";
		for (unsigned j = 0; j < 20; ++j)
			if (freq_matrix[j * 20 + j] != 0)
				os << valid_aa_chars[j]
					 << freq_matrix[j * 20 + j] << ",";
	}
}



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	lookup tables

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//________________________________________________________________________________________

//	initialize_freq_lookup_table

//________________________________________________________________________________________
void initialize_freq_lookup_table(short table[256][256])
{
	for (unsigned i = 0; i < 256; ++i)
		for (unsigned j = 0; j < 256; ++j)
			table [i][j] = 400;

	for (unsigned i = 0; i < 20; ++i)
		for (unsigned j = 0; j < 20; ++j)
			table [valid_aa_chars[i]][valid_aa_chars[j]] = i * 20 + j;

}

//________________________________________________________________________________________

//	initialize_bitscore_lookup_table

//________________________________________________________________________________________
void initialize_bitscore_lookup_table(char lookup_table[256][256], char gap_penalty[256][256])
{
	// bitscore lookup table
	memset (lookup_table, 0, 256 * 256 * sizeof (char));
	for (unsigned i = 0; i < sizeof (blosum_chars)-1; ++i)
		for (unsigned j = 0; j < sizeof (blosum_chars)-1; ++j)
		lookup_table[blosum_chars[i]][blosum_chars[j]] = BLOSUM_MATRIX[i][j];

	// gap extension penalty
	for (unsigned i = 0; i < sizeof(blosum_chars)-1; ++i)
	{
		lookup_table[blosum_chars[i]]['-'] = -1;
		lookup_table['-'][blosum_chars[i]] = -1;
	}

	// gap initiation penalty
	memset (gap_penalty, 0, 256 * 256 * sizeof (char));
	for (unsigned i = 0; i < sizeof(blosum_chars)-1; ++i)
		gap_penalty[blosum_chars[i]]['-'] = -11;
}

















//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	orthologs

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888




struct ortholog_region
{
	std::string				gene_id1;
	std::string				gene_id2;
	std::vector<t_segment>	loci_masked;
	std::string				description;
};
void get_gene_ids(std::vector<ortholog_region>& orthologs)
{
	std::string line;
	t_progress_indicator dots(cerr, 400);

	while (cin && getline(cin, line))
	{
		std::vector<std::string> data;
		// get maximum of four parameters
        perl_split (line, "\t", data, false, 4);
		if (data.size() < 4)
		{
			cerr << "\nData not four parameters: " << line << "[" << data.size() << "]\n";
			continue;
		}
		orthologs.push_back(ortholog_region());

		// save gene ids and description
		orthologs.back().gene_id1		= data[0];
		orthologs.back().gene_id2		= data[1];
		orthologs.back().description	= data[3];
		if (! data[2].length())
			continue;
		// split on commas and put into region
		std::vector<std::string> loci_str;
		perl_split(data[2],  "|",  loci_str);
		for (unsigned i = 0; i < loci_str.size(); ++i)
			orthologs.back().loci_masked.push_back(t_segment(loci_str[i]));


		++dots;
	}
}

bool expand_gaps(const string& gaps_str, string& prot, string& cdna, const string& gene_id)
{
	// empty string encodes no gaps!
	if (! gaps_str.length())
		return true;

	// split string of gaps
	vector<string> gaps;
	perl_split(gaps_str,  "|",  gaps);
	for (int i = gaps.size()-1; i >= 0; --i)
	{
		// gaps should not come after a protein or be humungously long
		unsigned beg = atoi(gaps[i].c_str());
		unsigned len = atoi(gaps[i].c_str() + gaps[i].find(':') + 1);
		if (beg > prot.length() || len > 80000)
		{
			return false;
		}
		prot.insert(beg, string(len,  '-'));
		cdna.insert(beg * 3, string(len * 3,  '-'));
	}
	return true;
}

//________________________________________________________________________________________

//	store_gene_ids_as_temp_table

//________________________________________________________________________________________
bool store_gene_ids_as_temp_table(std::vector<ortholog_region>& orthologs, t_pg_cursor& pg_db)
{
	if (!orthologs.size())
		return true;

	cerr << "Give to PostgreSQL the list of gene_ids whose sequence alignments are needed.\n";
	if (!pg_db.ExecCommandOk("CREATE TABLE t_gene_id(gene_id1 TEXT, gene_id2 TEXT, array_index OID)"))
	{
		cerr	<< "\nERROR!! There was an error CREATEing the table of gene_ids "
				<< "t_gene_id\n"
				<<	pg_db.ErrorMessage()
				<<	std::endl;
		return false;
	}

	t_pg_copy_guard copy_buffer(pg_db, "t_gene_id", cerr);
	if (!copy_buffer)
		return false;

	{
		// loop through each sequence
		t_progress_indicator dots(cerr, orthologs.size() / 50 + 1);

		for(unsigned i = 0; i < orthologs.size(); ++i, ++dots)
		{
			copy_buffer.buf += orthologs[i].gene_id1 + "\t" + orthologs[i].gene_id2 + "\t";
			ultoa(i, copy_buffer.buf) += "\n";
			copy_buffer.cached_write();
		}

		if (!copy_buffer.end_copy())
			return false;
	}

	if(!pg_db.ExecCommandOk("CREATE INDEX i_gene_id on t_gene_id(gene_id1)"))
	{
		cerr	<< "\nERROR!! There was an error CREATEing the index for gene_ids "
				<< "t_gene_id\n"
				<<	pg_db.ErrorMessage()
				<<	std::endl;
		return false;
	}

	return true;
}


struct t_error_dump_values
{
	const string& gene_id1;     const string& gene_id2;
	const string& prot1;        const string& prot2;
	const string& cdna1;        const string& cdna2;
	const string& gaps1;        const string& gaps2;
	unsigned align_beg;         unsigned align_end;
	t_error_dump_values(const string& gene_id1_,		const string& gene_id2_,
						const string& prot1_,           const string& prot2_,
						const string& cdna1_,           const string& cdna2_,
						const string& gaps1_,           const string& gaps2_,
						unsigned align_beg_,            unsigned align_end_):
						gene_id1(gene_id1_),            gene_id2(gene_id2_),
						prot1(prot1_),              	prot2(prot2_),
						cdna1(cdna1_),              	cdna2(cdna2_),
						gaps1(gaps1_),                  gaps2(gaps2_),
						align_beg(align_beg_),          align_end(align_end_){}

};
bool dump_params_to_errfile(ofstream& err_strm, const string& err_msg,
							const t_error_dump_values& dump)
{
	//length problem dump diagnostic
	err_strm<< "Error!!!\n\t"
			<< err_msg
			<< "This ortholog will be ignored:\n\t"
			<< "gene_id [ "
							<< dump.gene_id1		<< " / "
							<< dump.gene_id2		<< "]\n\t"
			<< "Gap strings: [ "
							<< dump.gaps1			<< " / "
							<< dump.gaps2			<< "]\n\t"
			<< "Protein lengths: [ "
							<< dump.prot1.length()	<< " / "
							<< dump.prot2.length()	<< "]\n\t"
			<< "Protein 1: [ " << dump.prot1		<< "]\n\t"
			<< "Protein 2: [ " << dump.prot2		<< "]\n\t"
			<< "cDNA lengths: [ "
							<< dump.cdna1.length()	<< " / "
							<< dump.cdna2.length()	<< "]\n\t"
			<< "cDNA 1: [ " << dump.cdna1			<< "]\n\t"
			<< "cDNA 2: [ " << dump.cdna2			<< "]\n\t"
			<< "Alignment [ "
							<< dump.align_beg		<< " / "
							<< dump.align_end		<< "]\n";
	return true;
}


//________________________________________________________________________________________

//	analyze_ortholog

//		die if return false

//________________________________________________________________________________________
bool analyze_ortholog(		ostream& os,
							std::vector<ortholog_region>& vec_orthologs, ofstream& err_strm,
							unsigned mask_region, short freq_lookup_table[256][256],
							t_pg_cursor& pg_db, unsigned row)
{
	string prot1 = pg_db.getvalue(row, 0);
	string prot2 = pg_db.getvalue(row, 1);
	string cdna1 = pg_db.getvalue(row, 2);
	string cdna2 = pg_db.getvalue(row, 3);
	string gaps1_str = pg_db.getvalue(row, 4);
	string gaps2_str = pg_db.getvalue(row, 5);
	unsigned align_beg = pg_db.getvalue<unsigned>(row, 6) - 1;
	unsigned align_end = pg_db.getvalue<unsigned>(row, 7) + align_beg;
	unsigned index = pg_db.getvalue<unsigned>(row, 8);

	if (index > vec_orthologs.size())
	{
		cerr	<< "\n\nError: The index returned [" << index << "]"
				<< " was greater than the rows of gene_ids [" << vec_orthologs.size() <<"]\n";
		return false;
	}
	ortholog_region& orthologs = vec_orthologs[index];

	// check gaps, alignment and lengths
	t_error_dump_values dump (orthologs.gene_id1,   orthologs.gene_id2,
							  prot1,                prot2,
							  cdna1,                cdna2,
							  gaps1_str,            gaps2_str,
							  align_beg,            align_end);
	if (!expand_gaps(gaps1_str, prot1, cdna1, orthologs.gene_id1) ||
		!expand_gaps(gaps2_str, prot2, cdna2, orthologs.gene_id2))
		return dump_params_to_errfile(err_strm, "Errors in Gap Lengths. ", dump);
	if (align_beg >= align_end || align_end > prot1.length())
		return dump_params_to_errfile(err_strm, "Alignment coordinates invalid. ", dump);
	if (prot1.length() != prot2.length() ||
		cdna1.length() != cdna2.length() ||
		cdna1.length() != prot1.length() * 3)
		return dump_params_to_errfile(err_strm, "Unmatched sequence strings. ", dump);

	//
	// merge overlapping regions
	//
	std::vector<t_loci> vec_loci(3);
	vec_loci[0].push_back(t_segment(align_beg, align_end));
	vec_loci[1].swap(orthologs.loci_masked);
	merge_overlapping_loci(vec_loci[1], vec_loci[2], align_beg, align_end);

	if (!(mask_region & eINVERSE))
		vec_loci.erase(vec_loci.begin() + 2);
	if (!(mask_region & eMASKED))
		vec_loci.erase(vec_loci.begin() + 1);
	if (!(mask_region & eWHOLE))
		vec_loci.erase(vec_loci.begin() + 0);
	assert(vec_loci.size());


	//____________________________________________________________________________________
	// get strings comprising concatenated domain/free regions
	//____________________________________________________________________________________
		vector<string> proteins;
		for (unsigned j = 0; j < vec_loci.size(); ++j)
			extract_seq_pair_with_loci(proteins, prot1, prot2, vec_loci[j]);

	//____________________________________________________________________________________
	// print gene_id
	//____________________________________________________________________________________
		os	<< orthologs.gene_id1		<< "\t"
			<< orthologs.gene_id2;


	//____________________________________________________________________________________
	// print protein lengths
	//____________________________________________________________________________________
		for (unsigned j = 0; j < proteins.size(); j += 2)
			if (!proteins[j].length())
				os << "\t\\N";
			else
				os << "\t" << proteins[j].length();

	//____________________________________________________________________________________
	// get protein percentage ids
	//____________________________________________________________________________________
		get_all_regions_percentage_id(os, proteins);



	//____________________________________________________________________________________
	// get Ka Ks
	//____________________________________________________________________________________
		vector<string> cDNAs;
		for (unsigned j = 0; j < vec_loci.size(); ++j)
			extract_cdna_pair_with_loci(cDNAs, cdna1, cdna2, vec_loci[j]);
		if (!get_all_regions_ka_ks(	os, err_strm, cDNAs))
			return false;


	//____________________________________________________________________________________
	// get nucleotide percentage ids
	//____________________________________________________________________________________
		get_all_regions_percentage_id(os, cDNAs);


	//____________________________________________________________________________________
	// get bitscore
	//____________________________________________________________________________________
		get_all_regions_bitscores(	os, prot1, prot2, vec_loci);

	//____________________________________________________________________________________
	// get mouse_to_human_matrix
	//____________________________________________________________________________________
		get_all_regions_matrix(	os, proteins, freq_lookup_table);

		os << "\t" << orthologs.description;

		os << "\n";
	return true;
}


//________________________________________________________________________________________

//	prepare_temp_directory

//		do yn00 in temp directory to avoid crud

//________________________________________________________________________________________
bool prepare_temp_directory(char* orig_directory)
{
	// save original directory for restore later
	if (!getcwd(orig_directory, 2047))
	{
		orig_directory[0] ='\0';
		cerr << "Error!\nCould not get currect directory.\n";
		return false;
	}

	// create temp directory if necessary
	struct stat temp_stat;
	bool dir_create = false;
	if (stat(tmp_file_path.c_str(), &temp_stat))
		dir_create = true;
	if (dir_create && ::mkdir( tmp_file_path.c_str(),
								S_IRWXU|S_IRWXG|S_IRWXO ) != 0 )
	{
		cerr << "Error!\n Could not create " + tmp_file_path + "\n";
		return false;
	}

	// change to temp directory
	if (chdir(tmp_file_path.c_str()))
	{
		cerr << "Error!\n Could not change to " << tmp_file_path;
		return false;
	}
	else
		cerr << "\tUsing " + tmp_file_path + " as the temp directory for yn00.\n";

	// write contro file for yn00
	ofstream f_yn00ctl("yn00.ctl");
	if (!f_yn00ctl)
	{
		cerr << "Error!\n Could not open " + tmp_file_path + "yn00.ctl\n";
		return false;
	}
	f_yn00ctl <<	"seqfile = seq.tmp\n"
					"outfile = yn\n"
					"verbose = 0\n"
					"noisy = 0\n"
					"icode = 0\n"
					"weighting = 1\n"
					"commonf3x4 = 0\n";
	f_yn00ctl.close();
	return true;
}

//________________________________________________________________________________________

//	analyse_all_orthologs

//		die if return false

//________________________________________________________________________________________
struct t_auto_restore_dir
{
	char dir_name [2048];
	t_auto_restore_dir(){dir_name[0] = '\0';}
	~t_auto_restore_dir(){	if (dir_name[0]) chdir(dir_name);}
};

bool analyse_all_orthologs(std::vector<ortholog_region>& orthologs, const string& sql_str,
							unsigned cnt_estimate,
							unsigned mask_region, short freq_lookup_table[256][256],
							const std::string& error_file,
							t_pg_cursor& pg_db)
{
	ofstream err_strm (error_file.c_str());
	if (!err_strm)
	{
		cerr << "\n\nError\n\tUnable to open " << my_path << "/lengths.err\n";
		return false;
	}


	// create binary cursor
	if (!pg_db.Declare(sql_str.c_str(), true))
	{
		cerr << "Failed to create a cursor to retrieve orthologuous sequences!\n";
		return false;
	}

	// save original directory for restore later
	t_auto_restore_dir save_directory;
	if(!prepare_temp_directory(save_directory.dir_name))
		return false;


	// time whole process
	t_portable_timer timer;
	{
		// get data in batches
		t_progress_indicator dots(cerr, 1);

		unsigned batch_size = (cnt_estimate / 50 + 1);
		while (pg_db.Fetch(batch_size))
		{
			int curr_batch_size = pg_db.Tuples();
			if (!curr_batch_size)
				break;
			dots++;

			// buffer output
			stringstream os;
			for (int i = 0; i < curr_batch_size; ++i)
				if (!analyze_ortholog(os, orthologs, err_strm, mask_region, freq_lookup_table,
										pg_db, i))
					return false;
			cout << os.rdbuf();
		}
	}
	double seconds = timer.reset();
	if (seconds > 120)
		cerr << "\tAnalysing orthologs took\t" << seconds / 60 << " minutes\n";
	else
		cerr << "\tAnalysing orthologs took\t" << seconds << " seconds\n";

	// close cursor
	if (!pg_db.Close())
		return pg_print_sql_err("close a cursor for retrieving ortholog sequences",
									pg_db, "CLOSE panda_cursor", cerr);
	return true;
}










//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	main


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main(int argc, char* argv[])
{
	// don't use C streams
	std::ios::sync_with_stdio(false);


	// get command line arguments

	// if no arguments print usage and exit
	if (argc == 1)
	{
		char * my_argv[2] ={ argv[0], "-h"};
		int argc = 2;
		t_program_args args;
		if (!process_command_line_options(argc, my_argv, args))
			return 1;
	}
	t_program_args args;
	if (!process_command_line_options(argc, argv, args))
		return 1;

	//____________________________________________________________________________________
	// 		Connect to PostgreSQL
	//____________________________________________________________________________________
	t_pg_cursor pg_db(args.db_login.c_str(), "panda_cursor");	//"dbname=test user=postgres"
	if ( pg_db.ConnectionBad() )
	{
		args.CERR().strm()	<< "\nERROR!! Connection to database ["
							<< args.db_login << "] failed.\n"
							<< pg_db.ErrorMessage() << "\n\n";
		return 1;
	}


	my_path = GetPathPart(argv[0]);
	paml_path = args.paml;
	tmp_file_path = args.tempdir.length() ? args.tempdir : my_path + "/yn00_tmp";
	yn00_exe_path = paml_path + " " + "yn00.ctl";
	cerr << "\tPAML control file from " << yn00_exe_path << "\n";



	short freq_lookup_table[256][256];
	initialize_freq_lookup_table(freq_lookup_table);
	initialize_bitscore_lookup_table(lookup_table, gap_penalty);

	//____________________________________________________________________________________
	// read ortholog data from STDIN
	//____________________________________________________________________________________
	std::vector<ortholog_region> orthologs;
	std::cerr << "Read supplied gene IDs...\n";
	get_gene_ids(orthologs);

	store_gene_ids_as_temp_table(orthologs, pg_db);
//	pg_db->Commit();
//	pg_db->BeginTransaction();

	//____________________________________________________________________________________
	// analyse each ortholog, retrieving the sequence from the specified database
	//____________________________________________________________________________________
	string sql_fields = "SELECT s1.sequence, s2.sequence, cdna1.cdna, cdna2.cdna, "
							"aln_gaps1, aln_gaps2, aln.bl2seq_align_start, "
							"aln.bl2seq_align_len, array_index ";
	string sql_tables = "FROM t_gene_id NATURAL JOIN orthologs." + args.taxon1 + '_' +
						args.taxon2 + "_alignment aln JOIN " +
						args.taxon1 + ".cdna_sequences cdna1 ON (aln.pi1 = cdna1.pi AND "
						"aln.acc_code1 = cdna1.acc_code) JOIN sequences s1 ON "
						"(aln.pi1 = s1.pi) JOIN " + args.taxon2 + ".cdna_sequences cdna2 "
						"ON (aln.pi2 = cdna2.pi AND aln.acc_code2 = cdna2.acc_code) JOIN "
						"sequences s2 ON (aln.pi2 = s2.pi)";

	if (!pg_db.ExecTuplesOk(("SELECT count (*) " + sql_tables).c_str()))
	{
		cerr << "Error:\n\tCould not retrieve the count of ortholog sequences to be "
					"aligned.\n\n";
		return 1;
	}
	unsigned cnt_estimate =	atoi(pg_db.getvalue(0, 0));
	cerr << "\t" << cnt_estimate << "\tortholog sequences to be processed...\n";


//	string sql_str = "SELECT sequence1, sequence2, cdna1, cdna2, "
//						"gaps1, gaps2, align_start, align_len, array_index FROM "
//						"t_gene_id NATURAL JOIN " +  args.table;
//	std::cerr << "Calculate ortholog statistics...\n";
//	if (!args.pg_db->ExecTuplesOk(("EXPLAIN  " + sql_str).c_str()))
//	{
//		cerr << "Error:\n\tCould not retrieve statistics from " << args.table << "\n\n";
//		return 1;
//	}
//	for (unsigned i =0, limit = args.pg_db->Tuples(); i < limit; ++i)
//		cerr << args.pg_db->getvalue(i, 0) << "\n";


	if (!analyse_all_orthologs(orthologs, sql_fields + sql_tables, cnt_estimate, args.region,
								freq_lookup_table, args.error_file,
								pg_db))
	{
		cerr << "Analysis failed ;-<\n";
		return 1;
	}


	cerr << "Success!\n";
	return 0;
}


